Identifying the impacts of extreme weather

Code

# Load libraries
library(terra)
library(sf)
library(tidyverse)
library(tmap)
library(raster)
library(stars)
library(janitor)
library(kableExtra)

Code

# Import night lights data
night_lights1 <- terra::rast(here::here("posts/2024-12-09-houston-blackout-visualization","data", "VNP46A1/VNP46A1.A2021038.h08v05.001.2021039064328.tif"))
night_lights2 <- terra::rast(here::here("posts/2024-12-09-houston-blackout-visualization", "data", "VNP46A1/VNP46A1.A2021038.h08v06.001.2021039064329.tif"))
night_lights3 <- terra::rast(here::here("posts/2024-12-09-houston-blackout-visualization", "data", "VNP46A1/VNP46A1.A2021047.h08v05.001.2021048091106.tif"))
night_lights4 <- terra::rast(here::here("posts/2024-12-09-houston-blackout-visualization", "data", "VNP46A1/VNP46A1.A2021047.h08v06.001.2021048091105.tif"))

Code

# Merge and process night lights data, cropping to the Houston area
houston_extent <- extent(c(-96.5, -94.5, 29, 30.5))

# Check if CRS matches before merging
if (crs(night_lights1) != crs(night_lights2)) {
  stop("CRS does not match between night_lights1 and night_lights2")
}
if (crs(night_lights3) != crs(night_lights4)) {
  stop("CRS does not match between night_lights3 and night_lights4")
}

night_lights_before <- terra::merge(night_lights1, night_lights2) %>%
  crop(houston_extent)
night_lights_after <- terra::merge(night_lights3, night_lights4) %>%
  crop(houston_extent)

# Create blackout mask
night_lights_change <- night_lights_after - night_lights_before
night_lights_change[night_lights_change > -200] <- NA


# Vectorize blackout mask
night_lights_change_poly <- as.polygons(night_lights_change) %>%
  st_as_sf() %>%
  st_make_valid() %>%
  st_transform(crs = "EPSG:3083")

Code

# Import infrastructure data
roads <- read_sf(here::here("posts/2024-12-09-houston-blackout-visualization",
                            "data", "gis_osm_roads_free_1.gpkg"), 
                 query = "SELECT * FROM gis_osm_roads_free_1 WHERE fclass='motorway'") %>%
  st_transform(crs = "EPSG:3083")

houses <- read_sf(here::here("posts/2024-12-09-houston-blackout-visualization",
                             "data", "gis_osm_buildings_a_free_1.gpkg"),
                  query = "SELECT * FROM gis_osm_buildings_a_free_1 
                         WHERE (type IS NULL AND name IS NULL)
                         OR type in ('residential', 'apartments', 'house', 'static_caravan', 'detached')") %>%
  st_transform(crs = "EPSG:3083")

Code

# Create road buffer and exclude highways from blackout mask
roads_buffer <- st_buffer(roads, dist = 200)
blackouts_200m <- st_difference(night_lights_change_poly, st_union(roads_buffer))

tm_shape(blackouts_200m) +
  tm_polygons() + 
  tm_layout(main.title = "Blackout Mask")

Create a set of maps comparing night light intensities before and after the storms

Code

# Calculate common min and max for both rasters after log1p transformation
min_val <- min(c(values(log1p(night_lights_before)), 
                 values(log1p(night_lights_after))), na.rm = TRUE)
max_val <- max(c(values(log1p(night_lights_before)), 
                 values(log1p(night_lights_after))), na.rm = TRUE)

par(mfrow = c(1,2))

# Plot before
terra::plot(log1p(night_lights_before),
     col = viridis::viridis(50),
     main = "Night Lights Intensity \n(February 7, 2021)\n",
     zlim = c(min_val, max_val),
    plg = list(title = "log1p(nW/cm²/sr)"))

# Plot after
terra::plot(log1p(night_lights_after),
     col = viridis::viridis(100),
     main = "Night Lights Intensity \n(February 16, 2021)\n",
     zlim = c(min_val, max_val),
     plg = list(title = "log1p(nW/cm²/sr)"))

Code

par(mfrow = c(1,1))

Create a map of the homes in Houston that lost power

Code

houses_filtered <- houses[blackouts_200m, ]

tm_shape(houses_filtered) +
  tm_dots(col = "#CBC3E3") +
  tm_compass(type = "8star",
             size = 0.8) +
  tm_scale_bar() +
  tm_layout(main.title = "Homes in Houston that Lost Power")

Estimate of the number of homes in Houston that lost power

Code

nrow(houses_filtered)

[1] 157970

Code

houses_count <- houses_filtered %>%
  group_by(type) %>%
  summarize(count = n()) %>%
  ungroup() %>%
  st_drop_geometry()

# Check that the sum of buildings adds up to the number of rows in `houses_filtered` 
testthat::expect_equal(sum(houses_count$count), nrow(houses_filtered))

houses_count_table <- houses_count %>%
  kbl(caption = "Number of homes that experienced power loss (total = 157967)") %>%
  kable_classic(html_font = "Cambria")

houses_count_table

Number of homes that experienced power loss (total = 157967)
type	count
apartments	1136
detached	353
house	19760
residential	1395
static_caravan	80
NA	135246

Create a map of the census tracts in Houston that lost power

Code

# Import census data
acs_texas <- st_read(here::here("posts/2024-12-09-houston-blackout-visualization",
                                "data", "ACS_2019_5YR_TRACT_48_TEXAS.gdb"), 
                     layer ="ACS_2019_5YR_TRACT_48_TEXAS")

Reading layer `ACS_2019_5YR_TRACT_48_TEXAS' from data source 
  `/Users/leilanie/git/leirubinstein.github.io/posts/2024-12-09-houston-blackout-visualization/data/ACS_2019_5YR_TRACT_48_TEXAS.gdb' 
  using driver `OpenFileGDB'
Simple feature collection with 5265 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -106.6456 ymin: 25.83716 xmax: -93.50804 ymax: 36.5007
Geodetic CRS:  NAD83

Code

acs_income <- st_read(here::here("posts/2024-12-09-houston-blackout-visualization",
                                 "data", "ACS_2019_5YR_TRACT_48_TEXAS.gdb"), 
                      layer = "X19_INCOME")

Reading layer `X19_INCOME' from data source 
  `/Users/leilanie/git/leirubinstein.github.io/posts/2024-12-09-houston-blackout-visualization/data/ACS_2019_5YR_TRACT_48_TEXAS.gdb' 
  using driver `OpenFileGDB'

Code

# Join the median household income from the previous 12 months to the census tract geometries
acs <- left_join(acs_income, acs_texas, join_by(GEOID == GEOID_Data)) %>%
  st_as_sf() %>%
  st_make_valid() %>%
  st_transform(crs = "EPSG:3083")

# Create a hounding box for the Houston area
houston_bbox <- st_bbox(c(xmin = -96.5, xmax = -94.5, 
                          ymin = 29, ymax = 30.5)) %>%
  st_as_sfc() %>%
  st_set_crs(4326) %>%
  st_transform(st_crs(acs))

# Filter to Houston area census tracts
houston_tracts <- acs[houston_bbox, ]

# Identify census tracts that that contained homes that experienced blackouts
houston_tracts_blackouts <- st_filter(houston_tracts, houses_filtered)

# Plot tracts
tm_shape(houston_tracts) +
  tm_polygons(col = "white") +
  tm_shape(houston_tracts_blackouts) +
  tm_polygons(col = "#CBC3E3") +
  tm_borders(col = "black") +
  tm_compass(type = "8star",
             size = 0.8) +
  tm_scale_bar() +
  tm_layout(main.title = "Census Tracts in Houston that Experienced Blackouts \n(shown in purple)",
            scale = 0.8)

Compare the distributions of median household income for census tracts that did and did not experience blackouts

Code

# Find census tracts that did NOT experience blackouts
houston_non_blackout <- houston_tracts %>%
  filter(!GEOID %in% houston_tracts_blackouts$GEOID)

# Label tracts by power status
houston_tracts_blackouts$status <- "Lost Power"
houston_non_blackout$status <- "Did Not Lose Power"
combined_tracts <- rbind(houston_tracts_blackouts, houston_non_blackout)

# Plot combined tracts
ggplot(combined_tracts) +
  geom_boxplot(aes(x = B19013e1, fill = status)) + # `B19013e1` = Estimate of median household income in the past 12 months (in 2010 inflation-adjusted dollars)
  labs(title = "Distribution of Median Household Income by Power Loss Status",
       x = "\nMedian Household Income (in 2011 inflation-adjusted dollars)",
       fill = "Power Status") +
  scale_fill_manual(values = c("Lost Power" = "#CBC3E3", 
                               "Did Not Lose Power" = "lightgreen")) +
  theme_bw() +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank())

Reflection

I found that median household income was slightly higher in the census tracts that lost power vs. those that did not lose power. Census tracts on the outskirts of the Houston area generally experienced less power loss, but most of the majority of tracts were still affected. According to Busby et al., all socio-economic groups were affected, but the effects of the storm were harder on low-income families, who live in older, poorly insulated homes and have limited access to resources for repairs and physical health after the storm.

Data Citations

Data	Citation	Link
Night Lights	NASA Earth Data (2024). Level-1 and Atmospheric Archive & Distribution System Distributed Active Archive Center (LAADS DAAC)	Level-1 and Atmospheric Archive & Distribution System Distributed Active Archive Center (LAADS DAAC).
Infrastructure Data	Planet OSM (2024). Retrieved from https://download.geofabrik.de/	Link
Census Data	U.S. Census Bureau. (2019). American Community Survey 1-year Public Use Microdata Samples. Retrieved from https://www.census.gov/programs-surveys/acs/news/data-releases.2019.html#list-tab-1133175109	[Link] (https://www.census.gov/programs-surveys/acs/news/data-releases.2019.html#list-tab-1133175109)
Understanding the 2021 Texas Blackout	Busby, J. W., Baker, K., Bazilian, M. D., Gilbert, A. Q., Grubert, E., Rai, V., Rhodes, J. D., Shidore, S., Smith, C. A., & Webber, M. E. (2021). Cascading risks: Understanding the 2021 winter blackout in Texas. Energy Research & Social Science, 77, 102106. https://doi.org/10.1016/j.erss.2021.102106	Cascading risks: Understanding the 2021 winter blackout in Texas